********************************************************************************
********************************************************************************
********* Tables for program-wide cost-effectiveness


********************************************************************************
********************** incorporating social cost of carbon

*******************************************************
***** table for program-wide CBA - sample with at least 12 months post
clear all
use "$maindir/Results/Model_Outputs/cba_scc.dta"
keep if nmonths==12

**** identifying the percentiles of net present benefits
sort year20_npv
cumul year20_npv, gen(pctile_year20)
replace pctile_year20 = round(pctile_year20, 0.01)
tostring pctile_year20, replace force format(%9.2f)

sort year30_npv
cumul year30_npv, gen(pctile_year30)
replace pctile_year30 = round(pctile_year30, 0.01)
tostring pctile_year30, replace force format(%9.2f)

sort year10_npv
cumul year10_npv, gen(pctile_year10)
replace pctile_year10 = round(pctile_year10, 0.01)
tostring pctile_year10, replace force format(%9.2f)

sort year40_npv
cumul year40_npv, gen(pctile_year40)
replace pctile_year40 = round(pctile_year40, 0.01)
tostring pctile_year40, replace force format(%9.2f)

sort disc6_npv
cumul disc6_npv, gen(pctile_disc6)
replace pctile_disc6 = round(pctile_disc6, 0.01)
tostring pctile_disc6, replace force format(%9.2f)

sort disc0_npv
cumul disc0_npv, gen(pctile_disc0)
replace pctile_disc0 = round(pctile_disc0, 0.01)
tostring pctile_disc0, replace force format(%9.2f)

foreach y in year20 year30 year10 year40 disc6 disc0 {
sort `y'_npv
gen obsid_`y' = _n
labmask obsid_`y', values(pctile_`y')
}


*** cumulative costs and benefits
foreach y in year20 year30 year10 year40 disc6 disc0 {
sort pctile_`y'
gen tot_costs = sum(Real_nonHS_cost)
count
sca nhomes = `r(N)'
egen tot_`y'_costs = max(tot_costs)
replace tot_`y'_costs = tot_`y'_costs/1000000
drop tot_costs

foreach x in 01 05 10 25 {
gen tot_costs = sum(Real_nonHS_cost) if pctile_`y'>="0.`x'"
count if pctile_`y'>="0.`x'"
sca nhomes`x' = `r(N)'
egen tot_`y'_costs`x' = max(tot_costs)
replace tot_`y'_costs`x' = tot_`y'_costs`x'/1000000
drop tot_costs
}

gen tot_benefits = sum(`y'_benefits)
egen cumul_`y'_benefits = max(tot_benefits)
replace cumul_`y'_benefits = cumul_`y'_benefits/1000000
drop tot_benefits

foreach x in 01 05 10 25 {
gen tot_benefits = sum(`y'_benefits) if pctile_`y'>="0.`x'"
egen cumul_`y'_benefits`x' = max(tot_benefits)
replace cumul_`y'_benefits`x' = cumul_`y'_benefits`x'/1000000
drop tot_benefits
}

gen cumul_`y'_npv = cumul_`y'_benefits - tot_`y'_costs 

foreach x in 01 05 10 25 {
gen cumul_`y'_npv`x' = cumul_`y'_benefits`x' - tot_`y'_costs`x'
}
}

** auxiliary variables to help determine cutoffs for program-wide cost-effectiveness
foreach y in year20 year30 year10 year40 disc6 disc0 {
gsort -`y'_npv
gen `y'_ce = sum(`y'_npv)
}


foreach y in year20 year30 year10 year40 disc6 disc0 {
	matrix input summary_mat = (.,.,.,.,.,.,.)
	foreach x in tot_`y'_costs cumul_`y'_benefits cumul_`y'_npv {
	sum `x'
	sca `x'all = `r(mean)'
	sum `x'01
	sca `x'01 = `r(mean)'
	sum `x'05
	sca `x'05 = `r(mean)'
	sum `x'10
	sca `x'10 = `r(mean)'
	sum `x'25
	sca `x'25 = `r(mean)'
	matrix mat`x' = (`x'all, `x'01, `x'05, `x'10, `x'25, . , .)
	matrix summary_mat = (summary_mat \ mat`x')
	}
	matrix summary_mat = summary_mat[2..., 1...]

	sort `y'_ce
	qui: sum obsid_`y' if `y'_ce>0 & `y'_ce[_n-1]<0
	sca numobs = `r(N)'
	if numobs==1 {
		sca obs_ce = `r(mean)'
		sca define pct_ce = 100 - real(pctile_`y'[`=scalar(obs_ce)'])*100
	}
	else if numobs==0 {
		sca define pct_ce = 100
	}
	local res = "`=scalar(pct_ce)'" + "\%"
	di "TB>TC at `res'"
	matrix summary_mat[3,6] = `=scalar(pct_ce)'

	gsort -`y'_npv
	qui: sum obsid_`y' if `y'_npv<=0 & `y'_npv[_n-1]>0
	sca obs_zeronet = `r(mean)'
	sca obs_zeronet = real(pctile_`y'[`=scalar(obs_zeronet)'])*100
	local res2 = "`=scalar(obs_zeronet)'" + "\%"
	di "MB=MC at `res2'"
	matrix summary_mat[3,7] = `=scalar(obs_zeronet)'

	matrix summary_mat = (summary_mat \ nhomes, nhomes01, nhomes05, nhomes10, nhomes25, . , .)

	matrix colnames summary_mat = "Full Sample" "Top 99% Homes" "Top 95% Homes" "Top 90% Homes" "Top 75% Homes" "TB>TC" "MB=MC"
	matrix rownames summary_mat = "Total Costs (million $)" "Total Benefits (million $)" "Net Benefits (million $)" "Number of Homes"

	esttab matrix(summary_mat, fmt("%9.2f %9.2f %9.2f %9.0fc")) ///
		using "$maindir/Results/Tables/cbatable_`y'_scc.tex", replace
}






********************************************************************************
********************** using retail energy prices

*******************************************************
***** table for program-wide CBA - sample with at least 12 months post
clear
use "$maindir/Results/Model_Outputs/cba_retail.dta"
keep if nmonths==12

**** identifying the percentiles of net present benefits
sort year20_npv
cumul year20_npv, gen(pctile_year20)
replace pctile_year20 = round(pctile_year20, 0.01)
tostring pctile_year20, replace force format(%9.2f)

sort year30_npv
cumul year30_npv, gen(pctile_year30)
replace pctile_year30 = round(pctile_year30, 0.01)
tostring pctile_year30, replace force format(%9.2f)

sort year10_npv
cumul year10_npv, gen(pctile_year10)
replace pctile_year10 = round(pctile_year10, 0.01)
tostring pctile_year10, replace force format(%9.2f)

sort year40_npv
cumul year40_npv, gen(pctile_year40)
replace pctile_year40 = round(pctile_year40, 0.01)
tostring pctile_year40, replace force format(%9.2f)

sort disc6_npv
cumul disc6_npv, gen(pctile_disc6)
replace pctile_disc6 = round(pctile_disc6, 0.01)
tostring pctile_disc6, replace force format(%9.2f)

sort disc0_npv
cumul disc0_npv, gen(pctile_disc0)
replace pctile_disc0 = round(pctile_disc0, 0.01)
tostring pctile_disc0, replace force format(%9.2f)

foreach y in year20 year30 year10 year40 disc6 disc0 {
sort `y'_npv
gen obsid_`y' = _n
labmask obsid_`y', values(pctile_`y')
}


*** cumulative costs and benefits
foreach y in year20 year30 year10 year40 disc6 disc0 {
sort pctile_`y'
gen tot_costs = sum(Real_nonHS_cost)
count
sca nhomes = `r(N)'
egen tot_`y'_costs = max(tot_costs)
replace tot_`y'_costs = tot_`y'_costs/1000000
drop tot_costs

foreach x in 01 05 10 25 {
gen tot_costs = sum(Real_nonHS_cost) if pctile_`y'>="0.`x'"
count if pctile_`y'>="0.`x'"
sca nhomes`x' = `r(N)'
egen tot_`y'_costs`x' = max(tot_costs)
replace tot_`y'_costs`x' = tot_`y'_costs`x'/1000000
drop tot_costs
}

gen tot_benefits = sum(`y'_benefits)
egen cumul_`y'_benefits = max(tot_benefits)
replace cumul_`y'_benefits = cumul_`y'_benefits/1000000
drop tot_benefits

foreach x in 01 05 10 25 {
gen tot_benefits = sum(`y'_benefits) if pctile_`y'>="0.`x'"
egen cumul_`y'_benefits`x' = max(tot_benefits)
replace cumul_`y'_benefits`x' = cumul_`y'_benefits`x'/1000000
drop tot_benefits
}

gen cumul_`y'_npv = cumul_`y'_benefits - tot_`y'_costs 

foreach x in 01 05 10 25 {
gen cumul_`y'_npv`x' = cumul_`y'_benefits`x' - tot_`y'_costs`x'
}
}

** auxiliary variables to help determine cutoffs for program-wide cost-effectiveness
foreach y in year20 year30 year10 year40 disc6 disc0 {
gsort -`y'_npv
gen `y'_ce = sum(`y'_npv)
}


foreach y in year20 year30 year10 year40 disc6 disc0 {
	matrix input summary_mat = (.,.,.,.,.,.,.)
	foreach x in tot_`y'_costs cumul_`y'_benefits cumul_`y'_npv {
	sum `x'
	sca `x'all = `r(mean)'
	sum `x'01
	sca `x'01 = `r(mean)'
	sum `x'05
	sca `x'05 = `r(mean)'
	sum `x'10
	sca `x'10 = `r(mean)'
	sum `x'25
	sca `x'25 = `r(mean)'
	matrix mat`x' = (`x'all, `x'01, `x'05, `x'10, `x'25, . , .)
	matrix summary_mat = (summary_mat \ mat`x')
	}
	matrix summary_mat = summary_mat[2..., 1...]

	sort `y'_ce
	qui: sum obsid_`y' if `y'_ce>0 & `y'_ce[_n-1]<0
	sca numobs = `r(N)'
	if numobs==1 {
		sca obs_ce = `r(mean)'
		sca define pct_ce = 100 - real(pctile_`y'[`=scalar(obs_ce)'])*100
	}
	else if numobs==0 {
		sca define pct_ce = 100
	}
	local res = "`=scalar(pct_ce)'" + "\%"
	di "TB>TC at `res'"
	matrix summary_mat[3,6] = `=scalar(pct_ce)'

	gsort -`y'_npv
	qui: sum obsid_`y' if `y'_npv<=0 & `y'_npv[_n-1]>0
	sca obs_zeronet = `r(mean)'
	sca obs_zeronet = real(pctile_`y'[`=scalar(obs_zeronet)'])*100
	local res2 = "`=scalar(obs_zeronet)'" + "\%"
	di "MB=MC at `res2'"
	matrix summary_mat[3,7] = `=scalar(obs_zeronet)'

	matrix summary_mat = (summary_mat \ nhomes, nhomes01, nhomes05, nhomes10, nhomes25, . , .)

	matrix colnames summary_mat = "Full Sample" "Top 99% Homes" "Top 95% Homes" "Top 90% Homes" "Top 75% Homes" "TB>TC" "MB=MC"
	matrix rownames summary_mat = "Total Costs (million $)" "Total Benefits (million $)" "Net Benefits (million $)" "Number of Homes"

	esttab matrix(summary_mat, fmt("%9.2f %9.2f %9.2f %9.0fc")) ///
		using "$maindir/Results/Tables/cbatable_`y'_retail.tex", replace
}
